When the harmonic trap potential is absent, i.e., γ = 0, we can get the following equilibrium points from Eqs. (5a)–(5d) by assuming R_x(0) = R_y(0) R₀, in which (ξ_x0,ξ_y0) and (p_x0, p_y0) are the initial position and initial momenta of the condensate: 6 ξ¯x=(lπ−py0)/ϕξ¯y=ξy0+(kπ+px0)/ϕ.

where l,k = 0,±1,±2,… When p_x0, p_y0, and ϕ are fixed, equation (6) gives multiple equilibrium points in (ξ_x,ξ_y) plane spaced by p/θ one by one in both ξ_x and ξ_y directions, which implies that there are multiply cyclotron center of the system. The more interesting is that the stability of the equilibrium points depends on l and k, which further determines the cyclotron dynamics of the condensate. Linear stability analysis shows that when l + k is even, the equilibrium points are stable, otherwise, when l + k is odd, the equilibrium points are unstable.^[33] Therefore, the distribution of the (un)stable equilibrium points in lattice plane changes periodically with the lattice site. So, the stability of the condensate will also change.

At this case, the analytical results of cyclotron orbit of the condensate can also be obtained: 7 cos⁡ϕ(ξx+py0/ϕ)+exp(ϕ2R028)cos⁡ϕ[ξy−(ξy0+px0/ϕ)]=−12exp⁡[(12R02+ϕ2R028)]H0MF,

where H₀^MF is the initial Hamiltonian of the system. We can find that the cyclotron orbit of the condensate given by Eq. (7) is a periodic function of ϕ. When the condensate is initially located at a stable equilibrium point, the condensate will be localized (see the first row in Fig. 1). When the condensate is initially located at an unstable equilibrium point, the condensate will diffuse disorderly with time (see the first row in Fig. 2). At this time, the stable orbit circling around the stable equilibrium point is only the single periodic orbit which is given by Eq. (7).

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. The dynamics of the condensate with ϕ = 0.126. From the first row to the second row: γ = 0,0.001, respectively. The cyclotron orbit (red lines) of the second row in the m–n plane is given by Eqs. (14) and (15). The initial conditions are (ξx0, ξy0) = (50,50),(px0 , py0) = (π,π).

When the harmonic trap potential is present, i.e., γ ≠ 0, the system (2) only has one stable equilibrium point, i.e., , and all the cyclotron orbits around this equilibrium point are stable. Assuming J = 1, then the following equations can be obtained from Eqs. (5c) and (5d): 8 p˙x=ξ¨x4e−1/R02−ξ˙x2, 9 p˙y=ξ¨y4exp[−(ϕ2R024+1R02)]−ξ˙y2−ϕξ¨x.

Combining Eqs. (5a) and (8) and Eqs. (5b) and (9), respectively, we have 10 ξ¨x=4e−1/R02−ξ˙x2[−ϕξ˙y−2γξx], 11 ξ¨y=4exp⁡[−(ϕ2R024+1R02)]−ξ˙y2[ϕξ˙x−2γξy].

Obviously, both the artificial magnetic field and harmonic confinement affect the cyclotron motion of the condensate. The coupling effects of the lattice potential, artificial magnetic field, and harmonic trap potential make the cyclotron dynamics of the condensate quite complexity. Although we cannot find the analytic result, we can obtain the asymptotical cyclotron orbit of the condensate under the condition of small cyclotron radius. By linearizing Eqs. (10) and (11) at the equilibrium point (0,0), we have 12 ξ∼¨x+4e−1/2R02γξ∼x+2e−1/2R02ϕξ∼˙y=0, 13 ξ∼¨y+4exp[−(12R02+ϕ2R028)]γξ∼y−2exp⁡[−(12R02+ϕ2R028)]ϕξ∼˙x=0.

Setting , the asymptotical orbit of the condensate center can be described by the following equations 14 ξ∼x=c12+c22cos⁡(ω1t−φ1)+c32+c42cos⁡(ω2t−φ2), 15 ξ∼y=(ω128γϕD1D2−ϕ2γ−12ϕD2)×ω1(c12+c22)sin⁡(ω1t−φ1)+(ω228γϕD1D2−ϕ2γ−12ϕD2)×ω2(c32+c42)sin⁡(ω2t−φ2),

where 16 ω1,2={2[(D1+D2)γ+D1D2ϕ2]∓2[(D1+D2)γ+D1D2ϕ2]2−4D1D2γ2}1/2, φ1=arctan⁡(c2c1),φ2=arctan⁡(c4c3), c1=16D1D22γϕsin⁡(py0+ϕξx0)+ω22ξx0[4D1(D2ϕ2+γ)−ω22](ω22−ω12)[4D1(D2ϕ2+γ)−ω12−ω22],c2=2{4D1D2γϕξy0−D1sin⁡px0[4D1(D2ϕ2+γ)−ω22]}ω1(ω22−ω12),c3=ξx0−16D1D22γϕsin⁡(py0+ϕξx0)+ω22ξx0[4D1(D2ϕ2+γ)−ω22](ω22−ω12)[4D1(D2ϕ2+γ)−ω12−ω22],c4=2{4D1D2γϕξy0−D1sin⁡px0[4D1(D2ϕ2+γ)−ω12]}ω2(ω12−ω22).

Here, parameters D₁ and D₂ express the effect of the optical lattice. When D₁ = D₂ = 1, the lattice is absent. We can see that the cyclotron dynamics of the condensate is affected not only by the artificial magnetic field ϕ and harmonic trapping γ but also by the lattice potential and initial conditions. On the other hand, the presence of the additional harmonic trap potential changes the equilibrium state properties and modifies the cyclotron dynamics of the system.

To understand the modification of the harmonic trap potential on equilibrium state and its stability of the system, figures 1–2 show the wave packet dynamics obtained by direct numerical simulations of Eq. (3) with (ξ_x0,ξ_y0) = (50,50), (p_x0, p_y0) = (π,π), R₀ = 5 for different ϕ. When ϕ = 0.126 and γ = 0, (ξ_x0,ξ_y0) = (50,50) is a stable equilibrium point. The condensate is localized at its initial position (see the first row in Fig. 1). At the same time, the amplitude of the condensate oscillates periodically with time. However, (ξ_x0,ξ_y0) = (50,50) with ϕ = 0.126 is not the equilibrium point when the harmonic trap potential is present. Due to the existence of the harmonic trap potential, i.e., γ = 0.001, there is only one equilibrium point of the system under the action of the optical lattices, artificial magnetic field, and harmonic trap potential. In this way, the condensate does cyclotron motion around the equilibrium point as time varies (see the second row in Fig. 1). The red lines in m–n plane are our asymptotic orbits given by Eqs. (14) and (15) and it is in good agreement with the direct numerical simulation of Eq. (3).

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. The dynamics of the condensate with ϕ = 0.189. From the first row to the second row: γ = 0,0.001, respectively. The cyclotron orbit (red lines) of the second row in the m–n plane is given by Eqs. (14) and (15). The initial conditions are (ξx0,ξy0) = (50,50),(px0, py0) = (π,π).

When ϕ = 0.189 and γ = 0, (ξ_x0,ξ_y0) = (50,50) is an unstable equilibrium point under the action of lattice potential. At this case, the cyclotron orbit is unstable and the condensate expand disorderly with time (see the first row in Fig. 2). But (ξ_x0,ξ_y0) = (50,50) with ϕ = 0.189 is not the equilibrium point when the harmonic trap potential is present. In this case, as shown in the second row in Fig. 2, the condensate does a stable cyclotron motion along the asymptotic orbit predicted by Eqs. (14) and (15). Therefore we can gain a conclusion that, the harmonic trap potential plays role in modifying the equilibrium state properties and stabilizing the cyclotron orbits of the condensate in optical lattice under the artificial magnetic field. In addition, the presence of the harmonic trap potential makes the orbit quite complexity, and the complexity of the orbits depends on the coupling of the lattice potential, artificial magnetic field, and harmonic trap potential.

In order to have a deeper understanding of the cyclotron dynamics of the condensate, the two eigenmodes of ω₁ and ω₂ given by Eq. (16) varying with ϕ are shown by the black lines in Fig. 3, where different lines represent different γ. When the harmonic trap potential is absent, i.e., γ = 0, ω₁ ≡ 0, the condensate does single periodic oscillations with frequency w₂/(2π). In this case, only the single periodic cyclotron orbit exists. However, when the harmonic trap potential is present, the cyclotron dynamics of the condensate is the composition of the two different eigenmodes, one is a low-frequency mode with frequency ω₁ → (2π), and another one is a high-frequency mode with frequency ω₂/(2π), which cause a series of complex dynamics phenomena. The frequency of the two modes increases with γ. As shown by the black lines in Fig. 3, the low-frequency mode ω₁ is almost independent of the lattice and decreases monotonously with ϕ, and ω₁ → 0 when ϕ → ∞. It means that for sufficiently larger artificial magnetic field (ϕ → π/2), the cyclotron dynamics of the system approaches to a single periodic orbit with eigenmode ω₂. However, when the lattice is present, the high-frequency mode ω₂ changes nonlinearly with ϕ, i.e., the high-frequency mode ω₂ drops down to the value of ϕ = 0 after it reaches the maximum, while when the lattice potential is absent, i.e., D₁ = D₂ = 1, ω₂ is monotonically increasing with ϕ (see the red lines in the Fig. 3). The lattice potential significantly modifies the mode w₂, which makes rich cyclotron dynamics of the condensate. We can get rich cyclotron orbits from Eqs. (14) and (15) which depend on the coupling effects of the lattice potential, artificial magnetic field, and harmonic trap potential. The cyclotron orbit of the condensate will be multi-periodic or quasi-periodic depending on ω₁/ω₂ being a rational number or not.

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. The frequency of the two eigenmode against ϕ for different γ. The different curves correspond to different γ. The black lines represent the case when the lattice potential is present. The red lines denote the situation when the lattice potential does not exist.

We can obtain single periodic orbit by setting c₁ = c₂ = 0 or c₃ = c₄ = 0 according to the asymptotical solution of Eqs. (14) and (15). Taking c₃ = c₄ = 0 as an example, γ and ϕ need to satisfy the relationship as follows: 17 D2ϕ2+γ−D2γD1+(γ+D2γD1+D2ϕ2)2−4D2γ2D1=2D2γϕξy0D1sin⁡px0, 18 D2ϕ2−γ+D2γD1−(γ+D2γD1+D2ϕ2)2−4D2γ2D1=2D2ϕsin⁡(py0+ϕξx0)ξx0.

So, the orbit of the condensate can be written as 19 ξx2+ξy2[−ω13/8γϕD1D2+ω1(ϕ/2γ+1/2D2ϕ)]2=c12+c22.

It shows that when γ and ϕ satisfy Eqs. (17) and (18), the trajectory of the mass center of the condensate is an ellipse. The cyclotron orbit is a single periodic orbit with frequency ω = ω₁/(2π) and period T = 2p/ω₁.

To confirm the above theoretical predictions, as an example, we set the initial position (ξ_x0,ξ_y0) = (−5,5) and the initial momenta (p_x0, p_y0) = (0.5,1.0). Then we can obtain the corresponding strength of the harmonic trapping γ = 0.018 and artificial magnetic field ϕ = 0.094 from Eqs. (17)–(18). The asymptotical orbit of the condensate given by Eq. (19) and the corresponding direct numerical simulation of Eq. (3) are presented in Fig. 4. Under such conditions, the condensate makes clockwise elliptical motion, because we choose p_x0 > 0 and p_y0 > 0. From Fig. 4 we can find that the asymptotical result is in good agreement with the numerical result.

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. The dynamics of the condensate and its contour plot with γ = 0.018, ϕ = 0.094. The elliptical orbit (red lines) is our theoretical prediction of Eq. (19). The initial position is (ξx0,ξy0 ) = (−5,5) and the initial momenta are (πx0, py0 ) = (0.5,1.0).

If ω₁/ω₂ is a rational number, the cyclotron trajectory given by Eqs. (14) and (15) will be multi-periodic. Setting ω₁/ω₂ = 1/Q, then from Eq. (16) we can obtain 20 γ=D2ϕ2[(Q2+1)/Q]D2/D1−1−D2/D1,

where 1/Q means a rational number. So given the initial conditions and the magnetic field intensity, we can get the corresponding strength of the harmonic trap potential which makes ω₁/ω₂ = 1/Q a rational number. We plot the orbits given by Eqs. (14) and (15) with ω₁/ω₂ = 1/2, 1/3, 1/4, 1/5 respectively in Fig. 5, where the initial position is (ξ_x0,ξ_y0) = (3,3). Using the strength of the artificial magnetic field ϕ = 0.15, the corresponding harmonic confinement γ = 0.0427,0.016,0.0095,0.0067 can be calculated from Eq. (20). The initial momenta are shown in Fig. 5. We predict that neutral atoms exposed to an artificial magnetic field and harmonic trap potential in optical lattices have various multi-periodic orbits as harmonic trap potential and initial momenta vary. Obviously, the orbit of the condensate is not only closed, but also has many intersections. That is to say, the orbit is multi-periodic. Under this condition, the larger the γ is, the less number of the intersection points is.

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. The cyclotron orbits in the xx–ξy plane with ϕ = 0.15 and (ξx0,ξy0 ) = (3,3). From the first row to the fourth row: γ = 0.0427,0.016,0.0095,0.0067, respectively. From the first column to the fifth column: (px0, py0 ) = (0.5,−0.5),(0.5,0),(0.5,0.5),(0,0.5),(−0.5,0.5), respectively.

We also plot the temporal evolution and the contour plot of the wave packet in Fig. 6 by direct numerical simulation of Eq. (3) for the cases of O1–O3 marked in Fig. 5. The red lines in m–n plane are our theoretical results. As time goes by, the mass center of the condensate is moving along the asymptotical orbit given by Eqs. (14) and (15) although the Gaussian shape of the condensate is no longer maintained. After one period, the mass center of the condensate returns to its initial position. The theoretical prediction is in good agreement with the numerical simulation.

Fig. 6.

	Figure Option ViewDownloadNew Window
	Fig. 6. The dynamics of the condensate and its contour plot given by numerical simulation of Eq. (3) for the cases marked by O1, O2, and O3 in Fig. 5. The red orbits are given by Eqs. (14) and (15).

If ω₁/ω₂ is an irrational number, the asymptotical cyclotron orbit of the condensate given by Eqs. (14) and (15) will be quasi-periodic. The presence of the lattice potential, artificial magnetic field, and harmonic confinement leads to abundant and complex dynamics phenomena. Figure 7 illustrates the asymptotical cyclotron orbits with (ξ_x0,ξ_y0) = (5,5) and γ = 0.0067 for different ϕ. In this case, ω₁/ω₂ is an irrational number. At the first column, the initial momenta (p_x0, p_y0) = (0,0), and at the second column, the initial momenta (p_x0, p_y0) = (π,π). In these cases, the cyclotron orbit of the mass center of the condensate is quite complex. The condensate is doing a non-repetitive cyclotron motion in an annular region. The trajectory of the condensate is not closed and intersected at some points, i.e., the orbit is quasi-periodic. We can find from Fig. 3 that ω₁ → 0 when ϕ = 0.98, so the cyclotron dynamics of the system approaches to a single periodic orbit with eigemode ω₂ (see the second row in Fig. 7). The temporal evolution of the condensate given by direct numerical simulation of Eq. (3) for the case of O4 marked in Fig. 7 is presented in Fig. 8, where red lines in m–n plane are our theoretical results. The numerical results are also agree with our theoretical results.

Fig. 7.

	Figure Option ViewDownloadNew Window
	Fig. 7. The cyclotron orbits in the xx–ξy plane with (ξx0,ξy0 ) = (5,5),g = 0.0067. In the first row ϕ = 0.252 and in the second row ϕ = 0.98. From the first column to the second column: (px0, py0 ) = (0,0),(π,π), respectively.

Fig. 8.

	Figure Option ViewDownloadNew Window
	Fig. 8. The dynamics of the condensate and its contour plot given by Eq. (3) for the case marked by O4 in Fig. 7.